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Abstract 

We derive exact expressions for the forces and torques between biaxial molecules 
interacting via the RE-squared potential, a recent variant of the Gay-Berne poten- 
tial. Moreover, efficient routines have been provided for rigid body MD simulations, 
resulting in 1.6 times speedup compared to the two-point finite difference approach. 
It has also been shown that the time cost of a MD simulation will be almost equal 
to a similar MC simulation, making use of the provided routines. 
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1 Introduction 



In molecular simulations, the van der Waals interactions have a prominent 
and essential contribution to the non-bonded interactions and are typically 
described using the Lennard-Jones(6-12) potential or its variants [1,2]. An 
interaction potential of this type between two extended molecules is assumed 
to be a double summation over the respective atomic interaction sites: 

U^nt{Ml,M2)= E E Ua{rij;t,j) (1) 

where Mi and A^2 denote the interacting molecules and Ua{-) is the atomic in- 
teraction potential, e.g the Lennard-Jones(6-12) potential. The required com- 
putation time for the exact evaluation of this double sum is quadratic in the 
number of interacting sites. In practice, a large distant interaction cutoff ac- 
companied by a proper tapering is used to reduce the computation cost. More 
sophisticated and efficient approximate summation methods such as Ewald 
summation and the Method of Lights [3] are also widely used. 

As an alternative approach. Gay and Berne [4] proposed a more complicated 
single-site interaction potential (in contrast to a more sophisticated summa- 
tion) for uniaxial rigid molecules which was generalized to dissimilar and bi- 
axial molecules by Berardi et al as well [5]. In response to the criticism of 
the unclear microscopic interpretation of the Gay-Berne potential [6], we have 
recently used results from colloid science to derive an interaction potential 
through a systematic approximation of the Hamaker integral [7] for mixtures 
of ellipsoids of arbitrary size and shape, namely the RE-squared potential. 
The parameter space of the RE-squared potential is almost identical to that 
of Berardi, Fava and Zannoni [5] , agrees significantly better with the numer- 
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ically evaluated continuum approximation of Eq. (1) and has no nonphysical 
large distant limit. It has been verified that the new potential is superior to 
the biaxial Gay-Berne potential in representing the atomistic interactions of 
small organic molecules as well [8]. Moreover, the potential of mean force is 
representable with the same functional form of the RE-squared potential with 
negligible error [8]. 

In an anisotropic coarse-grained potential model, a molecule (Ai) is treated 
and described like a rigid body, leading to a considerable speedup in numerical 
simulations while preserving the fundamental behavior of atomisic potentials. 
Neglecting the atomic details, each molecule is characterized by the position 
of its center (a vector r) and its orientation (a unitary operator A or a unit 
quaternion q). 

Due to the complexity of the functional form of such potentials, numerical 
finite differences are widely used for the evaluation of forces and torques 
in rigid body molecular dynamics simulations. The numerical differentiation 
methods are prone to round-off errors and are generally expensive in large 
scale simulations. 

In this article, we will derive analytic expressions for the forces and torques 
between two molecules interacting via the RE-squared potentials. A set of 
optimized routines will be suggested for an efficient implementation of the 
given expressions. Finally, a time comparison between the two-point finite 
difference and the analytic derivatives will be presented. 
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2 The RE-squared Potential 



As mentioned earlier, the RE-squared potential [9] is a coarse-grained descrip- 
tion of the attractive and repulsive interactions between two biaxial molecules. 
Each molecule is treated like a biaxial ellipsoid and is described by two char- 
acteristic diagonal tensors (in the principal basis of the molecule) S and E, 
representing the half radii of the molecule and the strength of the pole contact 
interactions, respectively. As mentioned earlier, the orientation of a molecule 
is described by a center displacement vector r and a unitary operator A, 
revolving the bases of lab frame to the principal frame of the molecule. 

The attractive and repulsive contributions of the RE-squared potential be- 
tween two molecules with a relative center displacement of ri2 = r2 — ri and 
respective orientation tensors Ai and A2 are respectively: 



^RE-squa.ed^^^^ A2, V,,) = (l + 3^712X12;^ 



X 

1 e=x,y,z \o'e^-\- /I12/2 



n n I ,) \ .J (2a) 



(Ai,A2,ri2) = ^(-j (l + -ryi2Xi2^jx 

n n : (2b) 

i=l e=x,y,z \ae + hu/bOa J 

where ^412 is the Hamaker constant (the energy scale) , Uc is the atomic interac- 
tion radius and a^^ and cr^*^ are the half-radii of ith ellipsoid (i=l,2). 7712 
and X12 are purely orientation dependant terms, describing the anisotropy of 
the molecules and hi2 is the the least contact distance between the ellipsoids. 

The structure tensor Sj and the relative well-depth tensor E, are diagonal in 
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the principal basis of ith molecule and are defined as: 

S, = diag{<7«,<j«,a«} (3a) 

E, = diag{E«,E»,E«} (3b) 

where E^^ and E^^ are dimensionless energy scales inversely propor- 
tional to the well-depths of the respective orthogonal configurations of the 
interacting molecules. For large molecules with uniform constructions, it has 
been shown [9] that the energy parameteres are approximately representable 
in terms of the local contact curvatures using the Derjaguin expansion [9,10]: 

E, = aediag(^,^,^| (4) 

The term Xi2 quantifies the strength of interaction with respect to the local 
atomic interaction strength of the molecules and is defined as: 

Xl2(Ai, A2, fi2) = 2f^2Br2^(Ai, A2)fi2 (5) 

where B12 is defined in terms of the orientation tensors Aj and relative well- 
depth tensors E^: 

Bi2(Ai, Ai) = A^EiAi + A^E2A2. (6) 

The term rji2 describes the effect of contact curvatures of the molecules in the 
strength of the interaction and is defined as: 

det[Si]/ai^ + det[S2]/ cr^ 

[det[Hi2]/((7i + (72)" 

Where ai is the projected radius of ith ellipsoid along ri2: 



7712 (Ai, A2, ri2) = -y^, {') 



a,(A„ ri2) = {r^,AjSfA,h2)-'^' (8) 
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and the tensor H12 is defined as: 

Hi2(Ai, A2, f 12) = - Af S?Ai + - A^S2A2. (9) 

(7i (72 

No trivial solution is available for the least contact distance between two 
arbitrary ellipsoids {hu) [6]. The Gay-Berne approximation [4,6] is usually 
employed due to its low complexity and acceptable performance: 

= I|ri2|| -f7i2, (10) 
where the anisotropic distance function ai2 [5] is defined as: 

^12 = Qrf2Gr2%2) ' (11) 
and the symmetric overlap tensor G12 is: 

Gi2 = AfS?Ai + Ai^S^A2 (12) 

We will also employ this approximation in our derivation and will omit the 

superscript GB for shorthand in the rest of the article. 

3 Analytic expressions for forces and torques 

The algebraic structure of the attractive and repulsive contributions of the RE- 
squared potential are essentially the same. Thus, both of the contributions are 
expressible with a proper template structure, defined as: 

t^r— - fi^r (1 + *.,.x.f ) X n n f^^yff-) 

V/2'12/ V riuy i=ie=x,y,z\ay +hi2/CaJ 

(13) 

We will work through this template in the derivation. One may yield to the 
explicit form of each of the contributions by giving appropriate values to the 
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a-superscripted parameters according to Eq. (2a) and (2b). 



In an interaction between the molecules A4i and M2, the exerted force and 
torque on the molecule is most easily evaluated by applying proper virtual 
displacements and infinitesimal rotations to the interaction potential. The 
exerted force and torque on Ali is trivially obtained using the third law of 
Newton, afterwards. 

We denote the first-order translational and rotational variation operators on 
the coordinates of A42 by 5t and Sr, respectively. The translational variation 
operator is formally defined on a scalar function F as: 



d 

St[F{Ai, A2, ri2); p, e] := e— F(Ai, A2, + ep) 



e=0 



(14) 



where the unit vector p points to the direction of translational variation and e is 
an (infinitesimal) scalar for which the variations are valued. The proper defini- 
tion of the rotational variation is more involved. The projection of the exerted 
torque on AI2 along the unit vector n is obtained by applying the infinitesi- 
mal orthogonal operator I + en.cr on the operator revolving the molecule from 
the lab frame to its current frame, i.e. A2 . The resulting orientation operator 
would be ^(I -I- en.cr) A^^ = A2 — eA2n.cr, according to the anti-symmetry of 
the principal rotation generators (cr). The discussion suggests the definition: 



5r[F{Ai, A2, ri2); n, e] := ex 

d 



— F(^Ai, A2 - eA2n,ri2 



e=0 



(15) 



where fl — h.a is the rotation generator corresponding to the direction n. 
Acting exclusively on the coordinates of the second molecule (^^2), the oper- 
ators may be used to define the exerted force and torque along p and around 
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n respectively as: 

SriUA + Un) 



Fm2,p = P (16a) 



Nx2,n = ^ — (16b) 



Applying either of the operators to the template potential Ua, we get: 

XTT ITT h ^12^X12 + ^^12X12 (n^ 

5Ua/Ua = CFha-r \ 7 ^"-12 1 



hl2 + cr6aXi2?7l2 V h 



il2 



2 



hi2 + har]i2Xi20- ^ § e=5,^ c^a^^ + hi2 ^ ^^^^ 
We will complete the derivation by providing explicit expressions for the first- 
order variations appearing in Eq. (17). 

3.1 Derivation of the first-order variations 

3.1.1 Rotational variation 0/7712 

Applying operator to 7712 (Eq. 7) and dropping off the constant terms, we 
arrive at: 



(^R^i2 1 Srg^ 1 5r det H12 5Ra2 2 det S2 



7712 2(Ti + (72 2 detHi2 (jf det Si/fTi + det 82/0-1 
The rotational variation of o"2 is: 



(18) 



^ = -\al5R (vl,AlS-^A2h2) (19) 

0-2 Z ^ ^ 

= ealrl^ ( A2fi)^S2 ' A2f 12 (20) 

We have used the symmetry and the anti-symmetry properties of and f2, 
respectively. The rotational variation of (5i?Hi2 is required prior to 5i?detHi2: 
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SrRu = -— (A^S^Asfl + (A2l2)^S^A2) - 

(72 ^ ^ 



^A-S^A, 

•^2 



It can be easily verified that, 



(21) 



dimr 



det[r + eA] = det[r] + e ^ detfr^^^] + 0{e^) 



(22) 



i=l 



for arbitrary T and A, where r*-*-' is defined as: 



kl 



Aj; k 



(23) 



An exphcit expression for SudetHu is feasible using Eq. (21) and Eq. (22). 
Equating F and A to H12 and dfiliu/e respectively, the first-order terms of 
Eq. (22) will evidently be equal to the variation we are looking for. 



3.1.2 Rotational variation of Xu 



Using Eq. (5), it is straightforward to show that: 



SrXi2 ^ rf2(5R(Bi2^)ri2 
X12 ffsB^a^fia 



(24) 



where: 



SrBi2 = -e 



(A20)^E2A2 + A^EaAafi 



(25) 



Using the mathematical relation: 



(r + eA)-^ = -er-^Ar-^ + Oie") 



(26) 



for infinitesimal e, together with Eq. (24) and (25) we finally reach to: 
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5rXi2 = SeffsB-' [(A2l2)^E2A2 + 



= 4eff2Br2' A^E2 A2f2Br2'ri2 



(27) 



where the symmetry of Ej, B12 and their inverses have been used. 



3.1.3 Rotational variation of hi2 

We will use the Gay-Berne approximation for the least constant distance de- 
fined by Eq. (10). Accordingly, the rotational variation of hi2 is a result of the 
change in the anisotropic distance function (T12: 



5Rhi2 = (^^(^12 - 0-12) = -5Rax2 



(28) 



where the rotationally constant term ri2 drops out. The term 5r(Ti2 is easily 
expressed in terms of 5rGi2'- 



5r'7\2 1 '5_Rri2G'i2^fl2 



(^12 



2 ri2Gi2^ri2 



(29) 



Eq. (12) together with Eq. (26) result in: 



5j^G^^ = eG^i [(A2n)^S^A2 + A^S2A2^^lGl-2^ 



-26Gr2' 



A^S^A20 



^12 



(30) 



where the symmetry of S2, G12 and their inverses have been used. Thus, we 
finally reach to: 



^Rhl2 — 2^'^12^12^12^-^2 S2A2S1G^2^^^12 



(31) 
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3.1.4 Translational variation ofrju 



The displacement of the molecule M.2 results in a change in the direction of 
the connecting vector r^. Up to first order, this change is expressed as: 



(0) , 

^(i)_ r\^'+^P 



I (0) , -II 

= ?S + —(p- {p.h2)h2) + 0{e') (32) 
ri2 ^ ' 

Defining a new auxiliary vector results in a cleaner derivation: 

u := P ' (33) 
r\2 

Accordingly, bvx2 is obviously e times u. Applying operator to 7712, we reach 
to: 



bTTl\2 1 ^l O"! + <)t<^2 1 dct H12 

(34) 



7712 2 (Ti + (72 2 det H12 

((^tCi) det Si/(Tj' + (5tO"2) det 82/(7! 



det Si/(7f + det S2/(7i 

We will follow the same steps as the rotational case. The translational variation 
of the projected diameter (7j is: 



= -\ea\ (u^ArSr2A,r,2 + fr2ArSr^A,u) 

= -e(7fu^AfSr^A,ri2 (35) 



It is also easy to verify that: 

(^rHi2 = -^A'f S?Ai - r^A'^ S^A2 (36) 



(7T 



Finally, we may express (5Tdet[Hi2] exphcitly using Eq. (22) in terms of H12 
and its translational variation, Eq. (36). 
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3.1.5 Translational variation of xu 
Applying 5t to Xi2, we get: 



Xi2 r?2Br2^ri2 

The numerator simplifies to: 



(37) 



5rri2Bi2^ri2 = e(^u^Bi2^ri2 + rf2Bi2^uj 

= 2eu^Br2'ri2 (38) 

We finally reach to: 

StXi2 = 4eu^Br2'ri2 (39) 
We have used the symmetry and the translational invariance of B12. 

3.1.6 Translational variation of hi2 

Both of the involving terms in the definition of hi2 contribute to (5t/ii2. The 
contribution of the center displacement is: 

(5^ri2 = eri2.p (40) 

and the variation of the anisotropic distance function may be expressed as: 

Sto-12 1 5Trl2Guri2 



U\2 2 v\2^Y^Vx2 

Expanding and simplifying the numerator, we reach to: 



(41) 



(5y(7i2 = --e(7^2U^ fl2 (42) 



Adding up the above contributions, we finally get: 

1 



S^hx2 = efi2.p + - 6(7^2 U G12 fi2 (43) 
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4 An Efficient Implementation for Rigid-body Molecular Dynam- 
ics Simulations 



Most of the required matrix and vectors products in the evaluation of the 
first derivatives using the given expressions will be already available once one 
gets through the evaluation of the interaction energy beforehand. Without a 
careful implementation, a minimal speedup is expected due to the consider- 
able redundancy of the algebra. Therefore, a proper integration between the 
variable spaces of all routines must be considered. The three provided rou- 
tines demonstrate a suggested implementation. The first routine evaluates the 
interaction energy while the second and third routines calculate the torque 
and force. The latter routines depend on portions of variable space of the first 
routine in order to skip the redundant matrix products. We have also omitted 
the e factors appearing in the variations beforehand as they will finally factor 
out, according to Eq. 16. In practice, one call of the first routine accompanied 
by three calls of each of the second and third routines are mandatory in order 
to evaluate the three components of the force and the torque vectors. We have 
compared the computation time of an efficient C-language implementation of 
the proposed routines [11] against a numerical two-point finite differentiation. 
A large scale comparison (Pentium-M 1.7Ghz, GCC4) indicates that the an 
evaluation of the interaction energy and the force and torque vectors takes 
38.6 IIS using the provided routines while the same calculation takes 62.2 /is 
with the finite difference approach, leading to 1.6 times speedup. Figures (1) 
and (2) have been drawn with the aid of the provided routines and show the 
typical behavior of interaction force and torque between two prolate molecules. 
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5 Monte-Ccirlo and Molecular Dynamics Simulations Time Cost 



Monte Carlo (MC) and Molecular Dynamics (MD) simulations are two main 
concerns in studying molecular systems. MC simulations are usually faster 
and more effective in the studying of steady states while MD simulations play 
a more prominent role in the studying of transition states of certain systems. 
Furthermore, there are certain cases where MC simulations are of little interest 
(specially where the dynamical behavior is demanded). 

A MC step is considered to be as revealing an a MD step once each molecule 
successfully move to a new position in the phase space. In a system of N 
molecules, each having an average number of M neighbors, the time consump- 
tion of a MC step roughly is: 

Tmc = aNM X te (44) 

where a is the inverse of the acceptance ratio (usually, a ~ 2 with a proper 
conditioning) and te is the average required time of an energy evaluation. The 
corresponding time consumption of an MD step would be: 

NM, , , , 

TmD = -^{-^E + TiT + Tt) (45) 

where tf and tt are the average excessive time required for a single force 
and torque evaluation in all three directions. Using the values obtained from 
a sample large-scale simulation (with an acceptance of 50%), the ratio of the 
time expenses turn out to be: 

Tmd_ ^ S8.6 if,s) ^ ^ ^ 

Tmc 2 X 2 X 8.9 (/xs) ' ^ ^ 

using analytical first derivatives. The same ratio would be 1.7 using finite 
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differences. Tfierefore, one will end up with a MD simulation almost as fast 
as a MC simulation using the provided analytical derivatives. 
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Routine 1: Evaluation of Ua 

1: ri2 <= ri2/||ri2|| 

2: Ti <= Af SfAi 

3: T2 <= A^S|A2 

4: s ^ (Fi + r2)-iri2 

5: (7i2 <S= l/-^|s^ri2 

6: Zi ■<= Aif 12 

7: Z2 A2ri2 

8: vi <J= Sr^Zi 

9: V2 <^= S^^Z2 
10: ai <= l/^^i 
11: (72 <^ l/^zfy^ 
12: H12 ^ ri/cTi + r2/cT2 
13: dn <^ det H12 
14: <s= det Si 
15: ds2 det S2 
16: A <^ dsi/al + ds2l(yl 
17: 1/ ^Jdn/iai + a2) 
18: w ^ (AfEiAi + A|^E2A2)-^ri2 
19: hi2 <^ ri2 - ai2 
20: r/12 <^ A/i/ 
21: X12 ^ 2rf2W 
22: Evaluate Ua using Eq. (13) 
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Routine 2: Evaluation of SnUa 
Require: Evaluated variable space of Routine (1) 
1: A <^ -A2(n.(7) 
2: p<{=Afi2 
3: 5^0-2 -crfp'^V2 

4: SrU^2 <= (Al^SiA + A^SiA2)/(72 - (5K(72/(7i)r2 

5: (^Rciff 4= 

6: for i = 1 to 3 do 

7: J <^ H^2^ {Defined in the corresponding section} 
8: SjidH <^ Sjidn + det J 
9: end for 

lU. dRr)i2 <= .2ia,+a2) ~ 2d^ A^f 

11: SrXi2 ^ -4w^A|'E2Aw 

12: Snhu <^ -|cr^2S^*s 

13: Evaluate SuUa using Eq. (17) 
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Routine 3: Evaluation of SrUa 
Require: Evaluated variable space of Routine (1) 
1: 7<j=p^ri2 
2: u<S= (p-7ri2)/||ri2|| 
3: Ui Aiu 
4: U2 <(= A2U 
5: 6t(Ti ^ ^CTiUi Vi 
6: 8t(J2 ^ — cr|u2 V2 



8: (^Tiii? 

9: for i = 1 to 3 do 

10: J <^ H12 {Defined in the corresponding section} 
11: ^xdn Srdn + det J 

12: end for 



7: StUu <= -^Fi 




13: StVu <^ r7i2 




2dH 




14: 5tXi2 4u^- 



w 



15: 5Thx2 <^ 7 + |(7i2u'^s 

16: Evaluate SrUa using Eq. (17) 
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Fig. 1. Typical interaction energy (red thick line) and torque (black thin line) 
between two prolate molecules (in arbitrary units), continuously evolved from 
side-by-side to cross configuration. The ellipsoids are identical, having half-radii 
[11:2: 0.5] (in arbitrary units) and vertically separated by 5 units of length. 
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Fig. 2. Interaction energy (red thick lines) and vertical force (black thin lines) be- 
tween two prolate molecules in two different configurations with respect to the 
vertical center separation (r). The ellipsoids are identical, having half-radii [11 : 
2 : 0.5] (in arbitrary units). The y-axis ticks correspond to the vertical force (in 
arbitrary units). 
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6 Conclusion 

We have derived analytical expressions for the forces and torques exerted on 
two molecules interacting via the RE-squared potential. Moreover, efficient 
routines have been provided for molecular dynamics simulations. A numerical 
investigation reveals that the provided routines are 1.6 times faster than a 
two-point finite difference approach. The evaluation of energy derivatives is 
the most expensive element in a MD simulation. Using the provided analytic 
derivatives, a MD simulation will run almost as fast as a similar MC simulation 
(Eq. 46). This speedup leads to the possibility of larger scale MD simulations of 
a wide range of materials such as liquid crystals and certain organic molecules. 
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